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Abstract. We study the effect of atomic relaxation on the structure of moire 
patterns in twisted graphene on graphite and double layer graphene by large scale 
atomistic simulations. The reconstructed structure can be described as a superlattice 
of ‘hot spots’ with vortex-like behaviour of in-plane atomic displacements and 
increasing out-of-plane displacements with decreasing angle. These lattice distortions 
affect both scalar and vector potential and the resulting electronic properties. At 
low misorientation angles (<^1°) the optimized structures deviate drastically from 
the sinusoidal modulation which is often assumed in calculations of the electronic 
properties. The proposed structure might be verified by scanning probe microscopy 
measurements. 
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1. Introduction 

The study of van der Waals heterostructures, artihcial materials made of stacked 
layers of different two dimensional (2D) materials, is a very promising new direction 
in nanoscience and nanotechnology [T]. Apart from perspectives in applications, such as 
vertical geometry tunneling transistors [2] , these systems are important for fundamental 
science. Graphene on hexagonal BN (h-BN), in particular, turns out to be an ideal 
playground to study statistical and quantum mechanics in tunable incommensurate 
potentials P 111 El El 0 [HI- When two lattices with different lattice constants or 
different orientations are stacked upon each other, a larger periodicity, known as a moire 
pattern, emerges. For graphene on h-BN, the moire pattern depends on both lattice 
constant mismatch and misorientation. It was recently shown both experimentally and 
theoretically that in this case an incommensurate-to-commensurate transition for small 
misorientations takes place PIHI. 

The case of identical lattices where misorientation is the only source of 
incommensurability may be physically different. The case of large misorientation 
30°) was studied intensively for graphite flakes on graphite after recognizing the 
drastic reduction in sliding friction often called superlubricity p. Conversely, stick-slip 
dynamics with high friction were observed for small misorientation angles (< 5°) and 
these angles were therefore less studied in this context. The case of small angles became 
interesting after experimental studies of misaligned double layer graphene m nn 
because of the strong effect on the electronic structure, such as the appearance of 
secundary Dirac cones and van Hove singularities [la Ea El ES]. In most theoretical 
works [la ES] the moire patterns were modeled by rigid misaligned layers and the 
effective potential on the electrons was considered as the superposition of crystalline 
potentials from the two rotated layers. 

In one dimension (ID), starting from the seminal papers by van der Merwe [T7] . 
many studies [18] have shown that the superposition of two slightly incommensurate 
periodicities lead to the appearance of a soliton lattice where commensurate (phase- 
locked) regions are separated by thin misfit regions. The 2D case is much less studied 
but a qualitatively similar behaviour, where vortices in 2D play the role solitons in ID, 
was found in models of adsorbed monolayers on surfaces ng. For the case of twisted 
double layer graphene, a density functional theory (DFT) calculation down to angles 
< 1° has been recently reported [20] • For the case of very small angles, a full DFT 
based minimization of the lattice was not performed in view of the diverging size of 
the supercell. It was assumed that the extrapolation to small angles of the calculated 
sinusoidal modulation for angles > 3°, could describe the out-of-plane distortions. 

Classical simulations for carbon have been shown to describe structural properties 
accurately [211 122] • They allow to study the very large samples needed for very small 
misalignment and can be used to verify the nature of the modulation in this limit. 
Here we use an atomistic model based on the REBO [21] and Kolmogorov-Crespi [2H] 
potentials to optimize the structure of graphene on graphite and of double layer graphene 
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as a function of the misalignment. 

We find an essential reconstruction of the geometrical moire pattern at very small 
angles. The modulation cannot be described by a sinusoidal function anymore as often 
assnmed PI- Instead, a 2D lattice of misfit regions (‘hot spots’) with large out-of-plane 
displacements separated by flat domains is formed. 


2. Model 


When two lattices with different lattice constants or different orientations are stacked 
upon each other, a larger periodicity emerges: a moire pattern. The lattice mismatch 
and relative orientation determine the size of the moire pattern. For the case of graphene 
on graphite, the lattice constants are equal and the size of the moire pattern can be 
written as PH EH] 


^lattice 

21 sin(d/2)| 


( 1 ) 


where 9 is the relative orientation of two layers, aiatuce = V^d with d the distance 
between carbon atoms. In figure [T]3,d is indicated by an arrow between the centers 
of the moire patterns. 

One may constrnct samples of twisted graphene satisfying periodic boundary 
conditions in the in-plane x, y directions by identifying a common periodicity PH ESI EH] 
in the two rotated layers as illustrated in figure [T^. For one layer, we define a supercell 
with a basis vector ti = nai + ma 2 where ai and 02 are the basis vectors of the graphene 
unit cell shown in fignre[pi with n and m integers and n > m> 1. For the second layer, 
a cell with same size and rotated by an angle 9 can be obtained by taking a basis vector 
t 2 = (n + m)ai — ma 2 . The common supercell is then obtained by rotating by 9/2 the 
cell with basis vector ti and by —9/2 the cell with basis vector t 2 - Each pair of n and 
m thus identifies a common supercell for two layers rotated by an angle 9 given by 

2n^ -|- 2nm — im? 

cos9 =—^ -(2) 

2{n^ + nm -\- m^) 

The supercell has sides of length 


^cell [f'll |f'2| ^lattice ^T Tim TfR 


(3) 


and contains N atoms, with 

N = + nm-\-mf). (4) 

For small angles, this procedure quickly results in very large supercells, see figure [TJd. 

The size and shape of the moire pattern does not change if the top layer shifts 
with respect to the bottom layer. In fact, a 60° rotation corresponds to a change from 
AA to AB stacking which can also be achieved by a translation. The symmetry of the 
hexagonal lattice combined with the translational freedom lead to a symmetry around 
30°. 
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Figure 1. (a) Schematic construction of the unit cell, (b) Number of atoms N for 
N <15000, n < 500 and m < n, mirrored around 30°. The dark blue line shows Nm 
and the dashed light blue line 3Nm- (c) Four supercells with one moire pattern in each 
of them ((n,m)=(7,6), Nm indicated by the filled dot in (b)) (d) A single supercell 
with three moire patterns in it ((n,m)=(17,l), Nm indicated by the empty dot in (b)). 
Note that a similar value of the distance between moire patterns a^, indicated by the 
arrow, can correspond to different supercells. 


The number of atoms in a single moire pattern is given by 

^ 


sin^0 


(5) 


This value also forms a lower bound for the number of atoms in the supercell since 
each supercell contains an integer number of moire patterns. The values of N on the 
minimum bound line in hgure[^ correspond to supercells with n — m= 1, like the one in 
hgure[T]3. The next highlighted line in figure corresponds to for supercells with 
m = 1 and three moire patterns inside them (see hgure [^) . The three moire patterns 
in the supercell are not exactly equal due to the discreteness of the lattice, in the same 
way as a shift of the top layer leads to a moire pattern that is not exactly the same as 
prior to the shift. In figure one can see that the center of one moire pattern is on an 
atom while the other is on the center of a hexagon. However, the differences between 
the moire patterns are tiny and cannot be distinguished experimentally m- 

The size of the moire pattern is thus determined by geometry, but its structure is 
determined by the atomic displacements that minimize the total energy. We have studied 
the effects of relaxation on the structure of the moire patterns of twisted graphene on 
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(a) (b) (c) 


Figure 2. The effects of relaxation are shown for a sample with (n,m) = (82,1), 
0 = 1.2° and a^=115.3 A. (a) The sample prior to relaxation, (b) the sample after 
relaxation, (c) The displacements of the atoms as the result of relaxation for a sample 
(n,m) = (17,1), 0 = 5.7° and a^=24.5 A. The colour indicates size and the arrow the 
direction of the displacements. 

a graphite substrate and of double layer graphene. The interactions between atoms 
within the layer are given by the REBO potential EB as implemented in the molecular 
dynamics code LAMMPS [23. The interlayer interactions are given by the registry- 
dependent Kolmogorov-Crespi potential without the normals [23], scaled to match the 
equilibrium distance d =1.3978 A of the intralayer potential. The combination of these 
potentials has been shown to accurately reproduce the potential energy surface of the 
graphite substrate [221, although the corrugation is possibly underestimated PI. 

Using the damped dynamics algorithm FIRE [SI] we relax the atomic positions to 
their minimum energy configuration. As a model of graphene on graphite we consider 
a mobile layer of graphene on a rigid substrate whereas to model double layer graphene 
we allow atomic relaxations in both layers. For graphene on graphite we consider not 
only full relaxation in all directions but also a case where the interlayer distance z is 
kept fixed at 3 A. As this value is smaller than the equilibrium value (3.32 A for AB 
stacking), the potential energy corrugation between layers is enhanced. 

3. Graphene on graphite 

We first examine the case where the distance z of the graphene layer to the rigid substrate 
is kept fixed and atoms are allowed to move only in the x and y directions. The A A areas 
become smaller while the AB stacked areas expand, as shown by the difference between 
figure]^ and figure [2)3. This effect is achieved by a small rotation of each moire pattern 
around the AA center as illustrated in figure where we show the displacements from 
the positions prior to relaxation. For clarity we show the displacement for a sample 
with ttjn = 24.5 A but similar vortex-like structures occur also for larger moire patterns. 

In figure we show the distribution of bond lengths for decreasing angles. The 
rotational pattern in the displacements can also be observed in the bond lengths: the 
triangular pattern for 9 = 2.1° and 9 = 1.2° changes to a windmill shape for the larger 
moire pattern with 9 = 0.46°. We see that by reducing the angle, the modulation of the 
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Figure 3. Bondlenghts of relaxed configurations for samples with z = 3 A.. The 
supercell is shown in black. The bottom panels show the bond length along the dashed 
diagonal line, (a) 9 = 2.1°, (n,m)=(47,l), =66.4 A. (b) 9 = 1.2°, (n,m)=(82,l), 

am =115.3 A. (c) 9 = 0.46°, (n,m)=(216,l), =302.6 A. 


bond length changes from an almost sinusoidal behaviour to rather localized contraction 
close to the AA areas. The ‘windmill’ patterns appearing at small angles resemble those 
found in uni as a function of the misfit energy, which in our case would correspond to 
an increasing penalty for AA stacking. 

The distribution of bond lengths observed at fixed z = 3 A is very similar to those 
found when atoms are free to move in z, although the differences are smaller than for 
z = 3 A, due to the increased distance and hence smaller corrugation. When the atoms 
are allowed to relax in the out-of-plane (z) direction, this is the main direction of the 
atomic displacements. In figure the average distance to the substrate z is shown, 
together with the minimum and maximum values. For 9 > 20°, the graphene lies fiat 
on the substrate. When 9 decreases, the minimum and maximum values of z are clearly 
different from the average. The spreading of z agrees with the DFT calculations of ra 
but the average value is not shown there. For 9 < 3°, the minimum and maximum 
value of z stay constant at the equilibrium values of the AA and AB stacking, but the 
average distance to the substrate decreases, implying that the AB stacked areas grow. 
The binding energy Eb = {2Egraphene — E)/N, shown in figure]^, stays constant for a 
large range of values as also found in Below 0 = 2°, a range not reported in ra. 
we find that the binding energy increases significantly. 

In figure we show the distance from the substrate z for three decreasing angles. 
In the panels showing the modulation along the cell diagonal, we see that the behaviour 
already found for the in-plane displacement is even more pronounced: for 9 = 2.1° 
the modulation is roughly sinusoidal whereas for 9 = 0.46° the ‘hot spots’ are very 
pronounced and separated by flat domains. 

The main effects produced by structural deformations on the electronic properties 
of graphene are the modulation of the pseudo-electrostatic potential V proportional to 
local uniform compression/dilatation and the appearance of a (pseudo)-vector potential 
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Figure 4. (a) Average height of the layer relative to the other layer after relaxation for 
various samples, (b) Binding energy per atom after relaxation. For a rigid AA-stacked 
sample at the equilibrium distance F^5=16.96 meV and for AB-stacked E}j=21.32 meV. 




Figure 5. Out-of-plane distance for samples where the graphene layer is relaxed in 
all dimensions. The bottom panels show the out-of-plane distance along the dashed 
diagonal line, (a) 0 = 2.1°, (n,m) = (47,l), =66.4 A. (b) 0 = 1.2°, (n,m)=(82,l), 

am =115.3 A. (c) 0 = 0.46°, (n,m) = (216,l), =302.6 A. 


A with components proportional to the shear deformations |32l|33]: 


Ax = ^ {ts - 


A, 


_ 1 

y — 2 + ^3 


^ 2 ) ~ (dt 

— 2tl) ~ —2l3tUxy 


( 6 ) 


where ~ 3 — 4 eV is the deformation potential for graphene, Ad = | (di + 0^2 + da) — d 
with di,d 2 ,d 3 the first neighbour interatomic distances in the deformed lattice, 
are hopping parameters in the deformed lattice, t ~ 3 eV is their unperturbed value, 
/? ~ 2 — 3 is the electronic Griineisen parameter. The resulting pseudo-magnetic field 
B = V X A, which has an opposite sign for the two valleys K and K’, can be estimated 

as |33ll3l] 


B 


he /2 N ^ 
e V3V dL 


(7) 


where u is a typical value of the shear deformation and L is a typical size of the spatial 



































Relaxation of moire patterns for slightly misaligned identical lattices 



(a) 0 = 5.7° (b) 0 = 0.46° 

Figure 6. u as in equationfor (a) 0 = 5.7° and (b) 0 = 0.46°. 


variation of the deformation which we can estimate as a^. 

To estimate the strength of the pseudo magnetic field, we calculate the magnitude 
of the shear deformations from the first neighbour interatomic distances in the deformed 
lattice as 

Y^3 (da — d2)‘^ + {d2 + ds — 2di)^ 

“=-2d-■ <*> 

Taking the parameters from figure |^,b we find in both cases B of the order of 
1 T, a value comparable to the one originating from ripples for graphene in Si02 m- 
The straightforward experimental evidence of the existence of this field would be the 
suppression of weak localization effects H. For the smallest angle the field is much 
less homogeneous. We have indicated the main directions of the vector potential in the 
areas of largest shear deformation by arrows. One can see that the pattern is not trivial 
and, for instance, cannot be described as superposition of the field created by magnetic 
fluxes. The amplitude of the modulation of V can be estimated to be a few meV. 

4. Double layer graphene 

In order to model double layer graphene, we no longer consider a rigid substrate, but 
two graphene layers which are both free to move in all directions. Whereas for graphene 
on graphite all corrugation occurs within one layer, in the case of a double layer the 
corrugation is shared between two layers, as shown in figure Furthermore we note 
that the deviation of the corrugation from the sinusoidal shape occurs for larger angles 
than in the case of graphene on graphite. Already at 9 = 1.2° the behaviour is no 
longer sinusoidal and the relative deformation of the two layers seems to lead to a more 
complex pattern than for a single layer. 

5. Conclusion 

We have shown that the lattice deformations in graphene on graphite and double 
layer graphene for small misorientation angle become very different from the sinusoidal 
modulation commonly used in theoretical models, leading to a structure of ‘hot spots’ 
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Figure 7. Out-of-plane distance for double layer graphene. The bottom four panels 
show ^ along the dashed line in the top figure. The dashed lines show the z for graphene 
on graphite as in figure 



separated by roughly uniform domains. The relatively large out-of-plane deformation at 
the Tot spots’ should be observable by scanning probe microscopy. Despite the modest 
effect of relaxation on the in-plane modulation, we estimate that it could give rise to 
a quite strong pseudomagnetic field, affecting the electronic structure of graphene. It 
would be interesting to pursue this topic further and study more precisely the effect of 
atomic relaxation of slightly misaligned graphene layers on the electronic properties. 
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